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ABSTRACT 

We have run high-resolution, three-dimensional, hydrodynamic simulations of the impact of comet 
Shoemaker-Levy 9 into the atmosphere of Jupiter. We find that the energy deposition profile is largely similar 
to the previous two-dimensional calculations of Mac Low & Zahnle, though perhaps somewhat broader in the 
range of height over which the energy is deposited. As with similar calculations for impacts into the Venusian 
atmosphere, there is considerable sensitivity in the results to small changes in the initial conditions, indicating 
dynamical chaos. We calculated the median depth of energy deposition (the height z at which 50% of the 
bolide's energy has been released) per run. The mean value among runs is « 70 km below the 1-bar level, for 
a 1-km diameter impactor of porous ice of density p = 0.6 g cm"^. The standard deviation among these runs is 
14 km. We find little evidence of a trend in these results with the resolution of the calculations (up to 57 cells 
across the impactor radius, or 8.8-m resolution), suggesting that resolutions as low as 16 grid cells across the 
radius of the bolide may yield good results for this particular quantity. 

Visualization of the bolide breakup shows that the ice impactors were shredded and/or compressed in a 
complicated manner but evidently did not fragment into separate, coherent masses, unlike calculations for 
basalt impactors. The processes that destroy the impactor take place at significantly shallower levels in the 
atmosphere (^ -40 km for a 1-km diameter bolide) but the shredded remains have enough inertia to carry them 
down another scale height or more before they lose their kinetic energy. 

Comparion basalt-impactor models shows that energy deposition curves for these objects have much less 
sensitivity to initial conditions than do ice impactors, which may reflect differences in the equation of state for 
the different kinds of objects, or a scale-dependent breakup phenomenology, with the preferred scale depending 
on impactor density. 

Models of impactors covering a ^ 600-fold range of mass (m) show that larger impactors descend slightly 
deeper than expected from scaling the intercepted atmospheric column mass by the impactor mass. Instead, the 
intercepted column mass scales as m'-^. 
Subject headings: hydrodynamics - comets: SL9 - planets and satellites: Jupiter 



L INTRODUCTION 

For the week starting 16 July 1994, fragments of comet 
Shoemaker-Levy 9 (SL9) hit Jupiter. Most of the world's tele- 
scopes observed the event, collecting an unprecedented vol- 
ume of imaging, photometry, and spectroscopy that spanned 
all sensible wav elengths. Papers in collections edited by 
IWest & Bohnhardt ( 1995) and Noll et al. ( 1996) review the 
data and early interpretation. Harrinaton et al. (2004) review 
the phenomenology, efforts to understand the impact phenom- 
ena, and open questions about the impacts. A brief summary 
of points relevant to this paper follows. 

The impacts followed the same basic phenomenology. The 
orbital path of the comet fragments intersected the planetary 
surface at ^^45° S latitude, and planetary rotation arranged 
for a girdle of well-separated impacts there. Each impactor 
fell into the atmosphere at over 60 km s"' and an impact 
angle of about 45° (Chodas & Yeoman s 1996). The ground 
track of the impactors moved toward the northwest. The 
bolides crushed, ablated, and decelerated as they fell through 
the atmosphere, leaving an entry channel filled with super- 
heated gas (e.g., Ahrens et al. 1994a b; Boslouah et al. 1995; 
Epslough & Gladstone 1997; Crawford 1996 ; Crawford et al. 
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'1994VMac Low' 1996^, 'Mac Low & Zahnle" 1994';'Takata et aj 
1994; Takata & Ahrens 1997; Shuvalov et al. 1999). This gas 
rushed back out the entry channel, exiting the atmosphere at 
an angle and flying ballistically into space. The Hubble Space 
Telescope (HST) resolved the plu mes, and saw several rise to 
--3000 km above the cloud tops ("Hammel et al."1995'). Ma- 
terial invisible from Earth rose higher still, as shown by its 
effects on the magnetosphere (Brecht et al. 2001, 1995). The 
plumes collapsed and re-entered the atmosphere in 20 min- 
utes, heating it and leading to infrared emission so strong that 
it was dubbed the "main event". In the initial hours after im- 
pact, two different syst ems of expandin g rings were seen, one 
in the visible by HST (Hammel et al. 1995) and the other in 
the infrared by McGregor et al. ( 1996). Peculiar patterns re- 
mained in the atmosphere, which were spread in latitude over 
the ensuing weeks ( Banfield et aljl99 CT. 

Chemistry was just as exciting. Main-event spectra were 
complex, but only five species were identified. Dozens ap- 
peared in the aftermath, and CO and CO2 are observed to this 
day, no w having crossed the equator into the no rthern hemi- 
sphere (iBezard et alJl2002t lLellouchet'ani2002h . Not to be 
outdone, the magnetosphere responded strongly and in a man- 
ner that changed throughout the week of impacts as ionized 
plume material loaded Jupiter's magnetic field lines. 

Initially viewed as perfect perturbation experiments, the im- 
pacts instead involved complex phenomena at fine spatial and 
temporal timescales. Models capable of reproducing many 
of the observed phenomena were too complex by several or- 
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ders of magnitude for the computational power available at the 
time. Many observers, including most spectroscopists, still 
have unpublished data that they cannot interpret. 

We made some headway in explaining phenomena 
of the images and lightcurves with simplified calcula- 
tions: Harrin gton & De ming (2001) extrapolated a pub- 
lished blowout velocity distribution to calculate the fluxes 
of mass, energy, and momentum on the atmosphere. 
[Deming & Harrington ( 2001 ) used those fluxes to drive a two- 
dimensional (2D) vertical-slice model of the atmosphere. The 
resulting lightcurves are a good match to data from 0.9 - 
12 [xm, and the time-dependent pressure (p) and temperature 
(T) grids show phenomena that mimic McGregor's Ring and 
other effects. 

Encouraged by this initial success, we have undertaken a 
study to model the impacts and their aftermath sufficiently to 
reproduce, in a realistic manner, all of the observed phenom- 
ena in the impact, blowout, plume flight, and plume splash 
phases. By "realistic", we mean that wherever we can elimi- 
nate an ad-hoc assumption or an analytic approximation, we 
do. We no longer use arbitrary initial conditions. Rather, 
we will use observations and the results of prior models for 
initialization, and will ultimately produce synthetic images, 
lightcurves, and spectra derived with radiative transfer from 
model results. For physical atmospheric effects, our ap- 
proach i s the chaining of hydrodynamic models outlined by 
iHarringt on et al. (2004). Our chemical models will be driven 
by tracer particle histories from the physical models. 

This paper presents the first results from our impact model. 
Since the observations of this phase were not as constraining 
as those of later phases, our primary goal was to produce a 
data grid with which to initialize subsequent models. Doing 
this believably required a look at, e.g., resolution convergence 
and sensitivity to initial conditions. Also, and as with any 
interesting investigation, there are serendipitous findings that 
re ach beyond this pa rticular set of i mpacts. 

iKorvcanskv et al.l (2002) and Ko rvcanskv & Zah nle (2003) 
have made 3D calculations of the impact of asteroids into 
the atmosphere of Venus using ZEUS-MP and its predeces- 
sor ZEUS3D. There are a number of similarities between the 
results of that work and the present study that we note below. 
Comparable hydrodynamical simulations of the initial phase 
of the impacts have bee n previously described by a num- 
ber of groups {Bosloug het all 119941: IC rawford et al. 1994; 



Grvaznov et al.' il994i iTak ata et al!" 1994; Yabe et al. 

Shoemaker et al. lll995t IC rawford et al. 1995; Svetsovl |l995[ 
^be et al. 1995; Crawford 1997; Shuvalov et al. 1997, 1999). 
Section 2 presents our model. We describe the results of 
over a dozen runs in section 3, and discuss their implications 
and future work in section 4. 

2. IMPACT MODEL 

For the calculations described in this paper we em ployed 
the ZEUS-MP hydrodynamics code ( Normanl EoOOh . which 
solves the equations for three-dimensional compressible gas 
flow. We have modified ZEUS-MP (base version 1.0) to 
include multiple materials. Modificat ions to t h e cod e are 
described in m ore detail ^by .Korvcanskv et alJ (120021) and 
IKorvcanskv & Z ahnle (20031). We have not included radia- 
tive transfer in ZEUS-MP; the short timescale and high opti- 
cal depth of the atmosphere below 1 bar, where the main dis- 
ruption of the impactors and energy deposition takes place, 
makes it unlikely that radiative transfer would significantly 
affect the dynamics. That assumption has been tested by 



IShuvalov et al.l ( Il999t . who indeed found that dynamical ef- 
fects of radiation transfer were insignificant, but that the im- 
pacting comets would be strongly heated at large depths. 

The Jovian atm ospheric profile comes from 
Deming & HarringtonI (|2p01). We also included minor 
modifications such as a moving grid that tracks the impactor 
at a variable velocity as it decelerates. Multiple materials 
were handled by the integration of tracer variables advected 
in the flow. For non-porous material, the tracer C gives 
the fraction of mass in the cell that is impactor material. 
Porosity is tracked by an additional tracer, and is treated with 
the so-called p- a m odel for a strengthless, porous solid 
dMenikoff & Kobeill999i) . The coordinate system (xi, X2, x^) 
is Cartesian and aligned with the bolide's initial velocity such 
that xi is the along-track coordinate, X2 is horizontal, and X3 
is perpendicular to the others. We relate (xi,X2,Xi) to local 
Cartesian coordinates {x,y,z) as follows: 

(1) 



X = X2 

y = -xi sinfl+jcscos^ 
z=xicos6+Xism6, 



(2) 

(3) 

where 9 is the angle between xi and the vertical. Note that 
local coordinates are not cardinal directions. Fluid velocities 
in the xi, X2, and xj, directions are v'l, V2, and V3, respectively. 
Other quantities that appear in the equations are the density, p 
and the internal energy per unit volume, e. The spatial resolu- 
tion of the calculations is described by the notation Rn, where 
n is the number of grid cells across the radius of the body 
in the high-resolution part of the grid. Away from an inner 
block of dimensions 4x2x2 km, the grid spacing increases 
geometrically (by a factor ^ 1 .04 per grid cell, depending on 
the overall resolution). The computational grid moves with 
the impactor and decelerates, keeping the object's front end 
about 1 km from the front end of the grid so that the object 
remains in the high-resolution portion as it disrupts. Calcu- 
lations in the paper were made with resolutions of R16, R32, 
and R57. 

We used the Tillotson equation of state, which was for- 
mulated for high- velocity impacts (Tillotson 1962; MelosQ 
119891) . although it cannot represent melting or mixed two- 
phase (gas-liquid) regimes. The Tillotson EOS has two 
regimes, one for cold and/or compressed material, the other 
for rarefied, hot conditions. For a given mass density of im- 
pactor material p and internal energy density (energy per unit 
volume) e, we have £ = e/p, rj = p/ Pt, ? = ??- 1, where pj is 
the density at zero pressure and £ is the specific energy. Then, 
for compressed states (77 > 1) and expanded ones (77 < 1) 
where £ < £ii,, the incipient vaporization energy, the pressure 
P is the sum of a thermal term and a cold-pressure term: 



P = p,= i a + 



Aq + Bq^ 



,e+— — (4) 

The parameters a, b, A, B, and £t are described in more detail 
and listed for common substances bv iMeloshI (11989): Kc is 
described below. For expanded states for which £ > £cv, the 
energy of complete vaporization. 



P = Pi,=ae + 



be 



Aq 



-5(I/r,-I)- 



(5) 



\+£/{£Tr]^) l+e-'^'i 

We modify the cold-pressure terms Aq + Bq^, and A^ in Eqs. 
^ and (|5} respectively, by a term of the form 1 + e~'^''' in the 
denominator, in order to provide a low-density cutoff as rec- 
ommended by .Melosh.(J982). We setKc = 1000. A strength- 
less material is assumed, so negative pressures (tensions) are 
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not sustainable; the pressure cutoff for low densities mimics 
that effect. For expanded states in which £iy < £ < £„., P is 
linearly interpolated between Pi and f/,. Then, on the grid, the 
resulting Tillotson pressure is weighted by the tracer C. The 
equation of state parameters for our materials (ice and basalt) 
are the same as used by Benz & Asphaug ( 1999) and are listed 
in table I. (As discussed below, we modeled basalt impactors 
in order to study the effect of the equation of state on impactor 
breakup.) 

For porosity, we use a simple model for a strengthless, 
porous solid based on the "p-a" formulation ( Her rmann! 
[1969; Menikoff & Kober 1999). We distinguish between the 
solid material density p„, (and energy e,„) and the same quan- 
tities for fluid in the porous void space {p„, and e,„). Space 
and material quantities are related by the porosity e, where 
0<e<l: p=(l- e)p„, and e = (1 - e)e„,. The material 
pressure p„ is provided by the equation of state PmiPm,em) 
(the Tillotson EOS described above) and the space pressure 
is thus p = (1 -e)p,„. Initially the porosity is l-(initial den- 
sity/nominal density) of the material. During the calculation, 
porosity is advected with the flow like the material tracer C. 
At each timestep the advected value of e is then adjusted by 
comparison with a function ef{p,pc) = eo(l -p/pcf" that de- 
pends on a "crushing pressure", pc, which for ice we set to 
10^ dyne cm"^. For p > pc, e/ = 0. The porosity is given 
by e = min(e,e/). During the calculation, porosity can only 
decrease from its initial value; if porosity is crushed out of 
a mass element, it does not return even if the pressure drops 
back below pc- 

The Tillotson EOS has been used for SL9 calculations by 
iGrvaznov et al. ( 1994) and Takata et al. ( 1994). Other equa- 
tions of state for cometary material used in previous calcu- 
lations include perfect gas (with adiabatic indices 7 = 1.2 
and 1.4), (Mac Low 1996), and a "stiffened gas" or Mur- 
naghan EOS with an additional perfect-gas thermal contribu- 
tion (Mac Low & Zahnle 1994). A similar stiff equation of 
state w as used by Yabe et al. (1994, 1995). Boslouah et al. 
( fT99l . ICrawford et al.. (.1994. J^95), and Crawford' (.1997.) 
used a tabular version of the ANEOS equation of state that 
took into account melting and vaporization. iShuvalov et al.l 
([1999) used a combination of a non-linear Gruneisen EOS 
in which the pressure increased roughly as the cube of the 
density ratio p/ pt combined with a tabular gas EOS for the 
vaporized component. 

The most important factors in the various equations of state 
are probably the relative stiffness dP /dp and the presence or 
absence of a dependence of pressure on the internal energy 
e. A relatively stiff EOS, and one in which P is also a func- 
tion of e, might be expected to result in calculations that show 
impactors blowing up at greater altitudes and exhibiting more 
radial spreading ("pancaking") than otherwise. This may ex- 
plain some of the differences among the results that have been 
found from different studies. We will not systematically ex- 
plore that question in this paper, but in view of the differences 
we found (discussed below) in results between ice and basalt 
impactors, we plan to address the issue in future work. 

Our "standard case" is a 1 -km-diameter, porous ice sphere 
impacting the atmosphere at v = 61 .46 km i"' and 9 = 43.09°. 
The velocity and imp act angle are the means of those of the 
21 comet fragments (IChodas & Yeomanslll996h in a frame 
rotating with Jupiter at the average latitude of the impacts 
(-44.02°). The gravitational acceleration at that latitude (in- 
cluding the J2 and centrifugal terms) is g = 2512 cm s~^. The 
bulk density of the bolide is p = 0.6 g cm"^ (lAsphaug & Benzi 



[T99i[T996HSolemll994l[T995l) . 

We have also carried out calculations of impacts with ob- 
jects of bulk density p = 0.1 and 0.92 g cm"^, the latter corre- 
sponding to non-porous solid ice. In addition, we have done 
calculations with impactors with volumes of 0.2, 3, 40, and 
125 times that of the standard case, or diameters of 0.584, 
1.44, 3.42, and 5 km. Not all combinations of density, diame- 
ters, and resolutions were run. 

One important result that emerged from the Venus- 
atmosphere calculations was the significant extent to 
which the results were sensitively dependent on initial 
conditions in a manner rem iniscent of dynamical chaos 
JKorvcanskv & Zahnlell2003h . Two calculations whose ini- 
tial conditions (for instance, impact velocity) differed by very 
small amounts (~ 0.1%) would develop large divergences in 
the course of their respective simulations. Korvcanskv et aD 
(|2000) took the view that the basic process of impactor dis- 
ruption was due to the growth of Rayleigh-Taylor and Kelvin- 
Helmholtz instabilities (Field & Ferrara 1995), and it is plau- 
sible that the seeds of the perturbations that grow to saturation 
are irregularities due to the finite resolution on the grid. In 
the physical case, one would expect analogous irregularities 
from the inevitable bumpiness of the bolide's surface. While 
the basic process (fragmentation and ablation) was always the 
same, there could be significant differences in the results of 
different runs for impactors of the same gross properties. For 
example, the diameter of the resulting crater on the surface of 
Venus might vary by several kilometers (on a scale of ~ 10 
km) depending on the exact details of how the event had un- 
folded. In the limit where the bolide just reached the surface, 
different cases might result in anything from a single crater 
several km in diameter, a group of small craters, or no crater 
at all. We might expect similar behavior to obtain in this case, 
and we have attempted to sample the distribution of results by 
performing several runs of the standard case with small dif- 
ferences in initial impact velocity Av of 0. 1 % of the initial 
velocity v (w 0.06 km s~') or displacements in the cross-track 
coordinates Ajc2, Ax3 by one half of a grid cell. For R16 cal- 
culations of the 1-km bolide the displacements are 15 m, for 
R32 they are 7.5 m, and for R57 they are 4.4 m. 

We use several diagnostics to characterize the calculations. 
The simplest and most significant is the profile of depo- 
sition of impactor kinetic energy dE/dz as a function of 
he ight z in the atmosphere, t he s ame quantity as discussed 
bv lMac Low & Zahnlel ( 11994 and'Mac Lowl ( 11996 *). The im- 
pactor kinetic energy E = E{z) was calculated by integrating 
the flux of kinetic energy passing through a given height z: 



E{z)= dt 



pC — ^- — V 1 dx2dx2 . 



(6) 



The area integral is taken over all cells whose height z = 
Xicos0 + X3sin0 is the desired value. [Note that the pro- 
jection factors cos6' in the element of area in the z plane 
(dA = dx2dxi/cos6) and the downward velocity (v^ = viCos6') 
cancel.] The time integral extends to the end of the calculation 
(typically 8-10 seconds). The density on the grid is weighted 
by the advected tracer field C that tags material that belongs 
to the impactor. We computed E{z) for -200 < z < 100 km 
at 1-km intervals. (For more massive bolides, E{z) was calcu- 
lated as far down as needed.) Having integrated E(z) for all 
heights, the deposition profile dE/dz follows by numerical 
differentiation. 

For a subset of the runs, we calculated spatially re- 
solved plots of mass flux, similar to those employed by 
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TABLE 1 
TillotsonEOS parameters 

material pr (gem"-*) a b A (erg c m~^) B (erg cm~^) Et (erg g~') g.yCerg g~') £,.,, (erg g~') 

JQg I 1 U I / I I ^ I I I U /I / ^ I ll 1 (' 

basalt 



TIT 



9.47 X 10'" 1.0 X 10" 7.73 X 10** 3.04 X 10 
2.67x10" 4.87 X 10'2 4.72 X 10'° 1.82 xlO" 



0.917 
2.70 



0.3 0.1 9.47 x 10' 
0.5 1.5 2.67 X 10" 



iJCorvcanskv & Zahnle ("2004") to study the fragmentation of 
asteroids in the Venusian atmosphere. Here we hope to corre- 
late events such as fragmentation, as revealed in the mass-flux 
plot, to specific features (like peaks) in the deposition profile. 
The time-integrated mass flux ^(z; jc,y) at a height z is given 
by 



^j.(z;x,y) = - I p{z;x,y)Ciz;x,y)vi(z;x,y)dt, 



(7) 



where z translates to the tilted plane Xicos0+X3sin0 in the 
grid coordinate system. In practice, /i is calculated as the 
accumulation of a set of integrals at time slices f,, each in- 
tegrated over the interval f; - Af,/2 to f; + Af,/2, by assuming 
that material moves at a constant velocity during that inter- 
val and counting all the mass that has crossed or will cross 
the plane during the interval. Due to the tilt of the grid, the y- 
location of the impactor is a function of z; in this case we refer 
the position of material on the plane to the centerline position 
defined by X2 = x^ = 0. 

The calculations were performed on a Beowulf cluster us- 
ing 32 2.4 GHz Opteron processors. The largest (R57) runs 
took ~ 3.1 X 10^ cpu seconds on a grid of 712 x 356 x 356 = 
9.0 X 10^ points. Impactors of 1 km diameter took '-^ 8- 10 
model seconds for the velocity to fall to 0.001 x the initial ve- 
locity, which was the condition for stopping the calculation. 
Timesteps were ^ 5 x 10"^ s during the initial phases of the 
impact, increasing to ^^ 3 x 10"^ s by the end, as the impactor 
slowed down. High-resolution output was written to disk ev- 
ery 0.1 model seconds, as noted above. The R57 calculations 
occupied 17 GB of memory and produced ~ 210 GB of out- 
put data. Not all data were saved from all runs. Wall-clock 
time for an R57 run was about three weeks. 

3. RESULTS 

Our main results are presented in Figs.[2and|2 which show 
the profile of kinetic energy deposition (dE /dz) for a number 
of realizations of the standard case described above. Figs. [2 
and|2]show much the same information, plotted in two differ- 
ent ways to emphasize two different points. Each panel shows 
profiles generated by calculations with very slight differences 
in initial conditions, as described above. The panels show cal- 
culations done at different resolutions: R16, R32, and R57. 
In Fig. [2 the horizontal scal e is logarithmic, to facilitate a 
comparison to th e results of .Mac Low & Zahnld tl994.) and 
lCrawfordl(IT996l) . 

For the most part previous calculations have been 2D ax- 
isymmetric calculations at moderately high resolutions (R25) 
with finite-difference (grid-based) hydrocodes. The exception 
is the calculation by Takata et al..Cl994.) . a 3D calculation em- 
ploying the smooth-particle hydrodynamics (SPH) method. 
Several groups found that a 1-km object penetrated more 
deeply than we found, to depths well below -100 km, P > 15 
bar (Boslouah et al. 1994; Crawford et al. 1994; Takata etaL| 
[1994; Crawford et al. 1995) . Shuvalov et al. ( 1999) found dif- 
ferent results, partly due to the different density and structure 



of the objects they model. Some of their calculations were 
made for an object of p=0.22 g cm"-'. They also investigated 
the impacts of objects of non-uniform density: a dense core 
(p = 1 g cm~^) and surrounding shell of low density (p = 0.1 
g cm""*), or a low-density object with high-density inclusions. 
Their objects had the same diameter (1 km) but were about 
1/3 as massive as our standard object. Shuvalov et al. found a 
rather broad, double-peaked distribution of energy deposition. 
Our re sults are most similar to those of Mac Low & Zahnla 
J 19941) and ICrawfordI Jl997h . The energy-release profile is 
sharply peaked, though not so strongly as found by Mac Low 
and Zahnle; a small amount of energy is deposited at levels 
below 100 km. It is apparent that the sensitivity to initial 
conditions described above also obtains for these simulations. 
This is brought out more strongly in Fig.|2] in which the hor- 
izontal scale is linear. 

An important question about simulations of this type is the 
degree of convergence that they exhibit as a function of nu- 
merical parameters such as resolution. In this case conver- 
gence means that extrapolation to infinite resolution of a series 
of models would yield a series of results that converged to the 
correct answer. Ideally, the limit of resolution-independent 
results is reached while we are still in the realizable limit of 
finite resolution. The question in this case is complicated by 
the sensitivity to initial conditions discussed above, since a 
degree of scatter in the results is introduced, as can be seen in 
Fig. 121 The scatter tends to obscure trends in the results as a 
function of resolution. We must thus compare the results as a 
function of resolution on a statistical basis. 

Also included among the runs plotted in Fig.|2]are two cal- 
culations (at resolutions R16 and R32) of bolides shaped like 
the asteroid 4769 Castalia but with the same volume and den- 
sity as our standard case. These runs are a test of the influence 
of the bolide's shape on the outcome, the object being in these 
cases an oblong object with axis ratios 2:1: 0.8 and perhaps a 
representative shape for non-spherical impactors. The object 
in these cases was oriented end-on to the direction of motion. 
The particular model shape was already available for use, 
having been extensively used in the calculations performed 
by Korvcanskv & Zahnle (2003). The results of the runs of 
Castalia-shaped impactors do not appear to be radically dif- 
ferent from spherical ones. Large changes in impactor shape 
do not seem to affect the outcome more than tiny changes 
in the initial position or velocity. Presumably, very extreme 
shapes (e.g., needle-like or plate-shaped objects) would have 
an effect, but moderately oblong shapes are not a strong influ- 
ence. For a relatively fragile impactor, it is probably true that 
the initial shape is quickly "forgotten" as the incoming object 
is rapidly deformed by aerodynamic forces, and that the same 
effects that apply to our spherical impactors also apply here. 
We expect that a sequence of calculations for irregular objects, 
with small changes in initial conditi ons, would show a simi- 
lar degree of scatter in the results. Crawf ord et alJ (11995 ) also 
simulated the impact of 2D, nonspherical objects (cylinders 
with length/diameter ratios of 1:3 and 3:1) and found a signif- 
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icant but not overwhelming dependence on body shape; para- 
doxically, the 3:1, rod-shaped impactor penetrated the least 
deeply among the three cases tested. 

Figure |3] shows various statistical measures of the energy 
deposition, in particular several different depths that charac- 
terize the results. Included in the plot are the mean depth 
(z = J z(dE /dz)dz/E) or the first moment of the energy de- 
position curve, the mode depth z, (the depth of peak energy 
deposition), and the depths, z„, at which n% of the bolide's 
energy has been deposited for n = 90, 50, and 10 (z5o is the 
median depth). The mean and median depths are quite simi- 
lar despite the skewness of the deposition profiles, while the 
peak energy-deposition depths z are somewhat deeper The 
trend (if any) of these measures as a function of resolution 
is weaker than the amount of scatter, as seen by the standard 
deviation. In particular the R57 runs give approximately the 
same results as the ones at lower resolution. Statistical tests 
(t- and F-tests) applied to the various measures of depth for 
R16 and R57 find that the differences between them are not 
statistically significant. However, the variance of the R16 re- 
sults is 3-4x larger than those of R57, and if the same results 
persisted after ^ 10 more cases were run, a significant result 
might emerge. That is, it is possible that the amount of scat- 
ter in the results is smaller for high resolution runs. We also 
note that the results for the Castalia runs do appear to result 
in systematically slightly shallower depths than the means at 
R16 and R32. However, only one Castalia run was done for 
R16 and R32, so the apparent results may not be significant. 

We also examined the spatially-resolved, integrated fluxes 
/i(z;x,y) described in the previous section. These pro- 
vide clues as to exactly how an impact proceeds, in terms 
of understanding the processes of fragmentation, ablation, 
and impactor spreading due to aerodynamic forces (pres- 
sure gradi ents). The last process has been denoted "pan- 
caking" (Zahnle 1992; Chvba et al. 1993) and has been mod- 
eled semi-analytically and applied to SL9 impacts by several 
groups (Chevalier & Sarazin 1994; Zahnle & Mac Low 1994; 
[ Mac Low & Zahnle. .1994: Jield &Ferrara.l995: .Crawford. 
Il997h . We will not make such a model in this paper, but sim- 
ply discuss the hydrodynamical results. 

Figure@]shows the time-integrated flux /i(z;x,y) calculated 
at various heights z relative to the 1-bar level for the R57 run. 
The flux is plotted on a logarithmic gray scale for /i > 5 x 10^ 
g cm"^, which ernphasizes relatively small values of /i. Values 
of 5 X 10"* g cm and above are shown in the deepest shade 
(black). (For comparison, note that a column 1 km high of 
p = 0.6 g cm"-' has a surface density of 6 x 10^ g cm"^.) Only 
the inner 4x4 km region of the projected grid is shown. Due 
to the zenith angle of the impact (« 44°) the "footprint" of the 
impact is elongated in the +y direction in the plots; the effect 
is most noticeable in the z = 20 km plot, in which the bolide is 
not yet strongly affected by the atmosphere. 

The object is torn apart quite high in the atmosphere (ap- 
proximately between z = -20 and -50 km) compared to 
the region below -50 km, where most of the energy is de- 
posited. Despite the explosion-like character of the impact, 
the crushed bolide has enough inertia to carry it down another 
scale height before it stops an d deposits its kinetic ene rgy. The 
same behavior was noted by IShuvalov et al.l Jl999h in their 
calculations. 

In Fig.|5]we show a similar quantity to /i for the R57 run, 
namely the column density a in the ixi,XT,) plane, or a "side 
view" of the impactor generated by integrating the density of 



impactor material in the X2 direction: 

a(xi,Xi) = / pCdx2. 



(8) 



Note that a is not time-integrated; we show particular instants 
in the calculation when the bolide is passing through z lev- 
els that are approximately the same as those shown in Fig.|3 
Note also that Fig. |5]is plotted in grid coordinates (xi, x^), 
which are related to z by Eq. |3] Figure |5] shows the com- 
pression and disruption of the bolide from a different point of 
view. The most interesting part of Fig.|5]is the initial compres- 
sion of the bolide seen at f = 2.7 s, followed by the shredding 
and sweeping back of material for f > 3 s. As noted above for 
Fig. @] the lateral spreading of the impactor during its initial 
compression is not larger than a factor of 2 or so. 

A notable feature is the character of the impactor dis- 
ruption. The impactor is apparently shredded and crushed, 
but does not seem to fragment into large pieces that sep- 
arate at significant velocity. Non-axisymmetric filamentary 
structures develop and then expand into a cloud of material. 
This behavior is different from what has been seen in cal- 
culations of asteroid impacts into the atmosphere o f Venus 
JKorvcanskv et al .'2002^.'Korvcanskv & Zahnle'2003') and in- 
ferred from craters on Venus ( Korvcanskv & Zahnle 2004) 
and the Earth ( Passev &M elosh 1980). Other calculations 
at lower resolution (R32 and R16) show similar behavior 
Given that the same code w as used for Venus calculations 
JKorvcanskv & Zahnlel2003l) and the calculations are in many 
ways very similar, we conclude that the material of the bolide 
(porous ice vs. rock) and its compressibility must control the 
character of impactor break-up in these physical situations. 

To test this idea, we ran an R32 simulation of a 1 -km spheri- 
cal impactor of non-porous basalt (p = 2.7 g cm"-') with other- 
wise identical conditions to our standard case. The results are 
shown in Fig.|6l where we compare an R32 porous-ice calcu- 
lation (top row) with the basalt impactor (bottom row) at se- 
lected heights in the atmosphere. Because the basalt bolide is 
~ 5 times as massive as the ice one, it penetrates more deeply, 
as reflected in the choice of heights in Fig. |6l Of more inter- 
est is the impactor breakup: the basalt object appears to break 
into several fragments and spread out (pancake) considerably 
more than does the ice impactor, whose degree of pancaking 
is rather modest, less than a factor of two in radius. However, 
to assess the degree to which the pancake model does or does 
not match the behavior seen in Figs.|3-|6li"equires quantitative 
modeling that we postpone to future work. 

Fig.Qshows an additional exploration of the possible out- 
comes due to differing substances. We ran five R16 impact 
calculations for 1-km diameter impactors of non-porous ice 
and porous and non-porous basalt. The differing impactor 
masses largely account for the variation of average penetra- 
tion depth per impactor type. However, the differences in 
the shapes and the level of variation among different trials 
is unexpected. The difference in outcomes between ice and 
basalt is marked; the icy impactors show much greater varia- 
tion in energy deposition (and a much greater overall spread 
in height) than do the basalt impactors. Differences in ini- 
tial porosity seem to have little effect compared with the dif- 
ference in material. Presumably this is due to differences in 
coefficients in the EOS. We speculate that the greater stiff- 
ness of the basalt EOS results in greater pancaking (as in Fig. 
|6j and (somewhat paradoxically) less variation in the disrup- 
tion and energy deposition as a result. One way to examine 
the question is to run models with simplified and artificially 
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Stiffened or softened equations of state. Understanding this re- 
sult would reveal fundamental impact physics and may enable 
simplified models of atmospheric impacts that do not require 
extensive hydrodynamic simulation. We hope to address this 
question in more detail in the future. 

As noted above, we have also run simulations of the im- 
pact of bolides of differing masses corresponding to masses 
0.2x, 3x, 40x, and 125x that of the standai'd case. The cor- 
responding diameters are 0.584, 1.44, 3.42, and 5 km. We 
ran 5 cases of each diameter at R16 resolution, perturbing the 
initial positions in X2 and xy, by one-half grid cell. The re- 
sults are shown in Fig.|8] where we plot the median depths of 
energy deposition zio, Zso, and zgo in the top panel as a func- 
tion of bolide mass. The bottom panel shows the same re- 
sult, but now we plot the corresponding atmosphere columns 
Aii0 50 90 = f P dz times the initial bolide cross-section 

A = 7rr^, normalized by the bolide mass m. Least-squares fits 
for z and fi are: 



Z50=-5.19 X 10-'(m/gf''**' km. 



n0.309 1 



Z90 

m 

m 
ficniA 



-1.82 X 10"V/g) km. 



= 8.16xlO-V/gy 
:1.20xlO-V/gy 
:1.03xlO"V/gy 



0.041 



0.199 



0.220 



(9) 



Zu) = -9.50 X 10-^(m/gy"'" km. 



Note that the same sensitivity to initial conditions appears to 
obtain across a 600-fold range in impactor mass. We expect 
that impactors would in general penetrate to column depths 
equivalent to their mass, i. e., nA ^ m , as was found roughly to 
be the case by Korvcanskv & Zahnld J2003h for impacts into 
the Venusian atmosphere. Thus, we would expect fiA oc m. 
As seen in Eq. |9] there is a weak dependence on impactor 
mass. More massive impactors penetrate somewhat more 
deeply than would be expected from a strictly proportional 
relation. 

4. CONCLUSIONS 
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To gain understanding of the SL9 impacts, we have carried 
out a number of 3D, hydrodynamic simulations of the impact 
of porous ice comets into the atmosphere of Jupiter We em- 
ployed the numerical hydrodynamics code ZEUS-MP, with 
some modifications to track the comet material (ice), its equa- 
tion of state, and the degree of porosity, if present. Calcula- 
tions were carried out at three different resolutions Rn, where 
n is the number of resolution elements across the bolide radius 
(for spherical impactors): R16, R32, and R57. We have paid 
special attention to the profile of energy deposition in the at- 
mosphere, as a measure of how deeply the bolide penetrated 
and for comparison with previous (mostly 2D) simulations. 
We carried out several calculations of a "standard case" (a 1- 
km-diameter, porous, ice comet with p = 0.6 g cm"-' and initial 
velocity like that of the SL9 impactors) with tiny variations in 
initial velocity or shifts in cross-wise position on the compu- 
tational grid, in order to test for sensitivity of the results to 
initial conditions, and to sample the distribution of results if 
the sensitivity were present. Such multiple calculations were 
carried out at all three of our resolutions to see if there was a 



convergence trend in the results. Two low and medium reso- 
lution calculations were also done of an impactor in the shape 
of the asteroid 4769 Castalia to see if there were noticeable 
differences for a non-spherical impactor. 

Energy-deposition profiles were fairly similar to 
those found for 2D calculations such as those done by 
Mac Low & Zahnle (1994), though they were slightly less 
sharply peaked. The median depth of energy deposition 
was « 70± 14 km below the 1-bar level, at an atmospheric 
pressure of «10 bars. The aforementioned sensitivity of 
the results to small changes in initial conditions produced 
significant variations in the energy deposition profiles, and 
the error bar just given refers to standard deviation of the 
individual profiles. Comparing the results of calculations at 
different resolutions shows very little trend in the results, 
compared to their scatter This suggests that, for the pur- 
pose of determining the depth of penetration of impactors, 
relatively low resolution (R16) is sufficient. 

We have visualized some of the calculations to learn about 
the impact process and how a bolide responds to aerodynamic 
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forces. Th e pictures we s ee ar e consistent with th e "pancake 
model" of Zahnle ( 1992) and Chyba et al. ( 1993). The im- 
pactor is flattened quite strongly at early times, but the extent 
of radial spreading was no more than a factor of two in radius. 
Shortly thereafter, the impactor develops non-axisymmetric 



structures and shreds into filamentary structures before com- 
ing apart completely. Material is blown back and outward 
as ablation proceeds until the impactor material expands into 
a cloud that slows down and deposits its kinetic energy into 
the atmosphere. The disruption takes place at a considerably 
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shallower depth (at ^ -40 km) than the peak deposition of en- 
ergy; the broken-up impactor material has sufficient inertia to 
be carried downward a significant distance (^ 1 scale height 
or more) before being stopped. 

A set of "low-resolution" (R16) runs of impactors over a 
600-fold range in mass (corresponding to diameters 0.584 < 
d <5 km) produced median depths of energy deposition rang- 
ing from 35 km to ^ 250 km below the 1-bar level. Scaling 
the results by the amount of atmospheric mass intercepted by 



the bolide showed a weak dependence on impactor mass, with 
more massive bolides penetrating slightly more deeply than 
predicted by a linear relation between intercepted atmospheric 
column and bolide mass. 

Future work will extend these results in a number of direc- 
tions. We will explore the parameter space of impactor mass 
and composition. High-resolution calculations will also serve 
as the basis for new models of the impactor plume and splash- 
back that generated the greater part of the SL9 phenomena ob- 
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served from Earth. The results of these calculations will also 
serve as input for simplified, semi-analytic, general models 
of atmospheric impacts for diverse situations such as impacts 
into the atmospheres of Earth, Venus, and Titan. 

We thank K. Zahnle for useful discussions, and the referee. 
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